Space representation of stochastic processes with delay 
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We show that a time series xt evolving by a non-local update rule xt = f(xt- n ,Xt-k) with two 
different delays k < n can be mapped onto a local process in two dimensions with special time- 
delayed boundary conditions provided that n and k are coprime. For certain stochastic update rules 
exhibiting a non-equilibrium phase transition this mapping implies that the critical behavior does 
not depend on the short delay k. In these cases, the autocorrelation function of the time series is 
related to the critical properties of directed percolation. 



I. INTRODUCTION 



Dynamical systems with time-delayed feedback show 
interesting phenomena which are relevant to a broad 
spectrum of research fields, such as nonlinear dynam- 
ics fiL laser physics [2j, neurobiology [3], chaos con- 
trol [4], synchronization jBI, @] and communication 0]. 
The mathematics of delayed differential equations is less 
developed and not as understood as that of ordinary 
ones [8]. Hence the theoretical analysis of delayed feed- 
back currently attracts a lot of attention Q. 

Most of the research on time delayed systems concen- 
trates on deterministic systems. However stochastic sys- 
tems with time-delayed feedback are discussed as well, 
for instance, in the context of gene regulation [101 ]. 

In this paper we investigate a simple discrete model 
with delay: a stochastic process for a single binary vari- 
able which evolves according to its own history. We show 
that this model can be mapped onto a two-dimensional 
stochastic cellular automaton in such a way that the 
time-delayed couplings become local. As a result, the 
autocorrelation function of the corresponding time series 
is related to the critical properties of directed percola- 
tion. Our results are however more general and may be 
applied to a large variety of systems. 



II. REORDERING OF THE TIME SERIES 



We consider a time series x t with a discrete time vari- 
able t G N that evolves by non-local stochastic updates 
(x t -n->%t-k) —> %t in such a way that the probability for 
the outcome x t is given by 



P(x t ) = f(x t ,Xt-n,X t -k) , 



t> n> k 



(1) 



where n and k are two different delays. The initial con- 
figuration may be given by specifying n subsequent ele- 
ments of the time series, e.g. xo, . . . , x n -±. The type of 
data represented by x t as well as the function / is not 
restricted in any way; the only important ingredient is 



that a new entry of the time series x t depends on the 
previous values at times t — n and t — k. 

Here we show that irrespective of the structure of x t 
and /, it is possible to rearrange the time series in such a 
way that the couplings become local in a two-dimensional 
representation. More specifically, we show that it is possi- 
ble to define a reordered series y t that evolves by updates 
(yt-niVt-s) -> Vt with 



p (yt) = f(yt,yt- n , yt- 3 ) , 



(2) 



where s is a different delay that switches between the 
values 1 or n + 1. The only condition for this transfor- 
mation to work is that the original delays n and k have 
to be coprime, i.e., they have no common divisor other 
than 1. The reordered time series y t is related to the 
original one by 



Xt = 2/T(t) 

where the map T(t) is given by 

T(t) = (tq) mod n + n[t\ n . 



(3) 



(4) 



In this equation [t\ n denotes integer division by n while 
q G {0, . . . , n — 1} is an integer such that 



(kq) mod n = 1 



(5) 



Note that the existence of q requires n and k to be co- 
prime. The condition of coprimality of n and k guaran- 
tees that the mapping is one-to-one with the inverse 



T _1 (Y) = (tk) mod n + n[t\ n . 



(6) 



Coprimality should not be seen as a shortcoming of the 
mapping since the number of coprimes to a given n, as 
given by the Euler's totient function 0(n), is known to 
increase sufficiently rapidly (faster than y/n) so that in 
the large n limit the transformation can be applied to 
systems with various delays k. 

To understand this transformation, the time line 
should be thought of as being divided into equidistant 
blocks of size n. The transformation simply reorders the 
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FIG. 1: Updating scheme in systems with delays n — 7 and 
/c = 3 before (upper figure) and after (lower figure) the trans- 
formation (|4j). The arrows indicate the couplings needed to 
update the variable X7 as a function of xo and X4- 

time indices in each block. This reordering takes place 
in the first term on the r.h.s. of Eq. (|4j) while the second 
term simply enumerates the blocks in such a way that 
they do not mix. 

As an example let us consider the special case n = 7 
and k = 3, hence q = 5. One starts out with a given 
configuration {xo, xi, . . . ,xq}. From this initial condition 
the whole time series can be constructed by iteration of 
Eq. (pQ). Applying the transformation Eq. (jU the first n 
sites are reordered by 

0^0; 1^3; 2^6; 3^2; 4^5; 5^1; 6^4. 

(7) 

The same reordering scheme takes place in the subse- 
quent blocks. As illustrated in Fig. HJ this transforma- 
tion preserves the long delay n while the short delay k is 
mapped onto 1 or n + 1. In the appendix we prove that 
this is also true in the general case as long as n and k are 
coprime. 




FIG. 2: (Color online) Original time series Xt with delays 
n — 7 and k = 3 in a two-dimensional representation [ill ]. 
The time series is divided into segments of size n and plot- 
ted line by line on top of each other. Starting with a given 
configuration {xo, xi, . . . , xq} the whole time series can be 
constructed by iteration of Eq. JT]), where the black arrows 
represent the long delay between and x t while the red 

ones indicate the short delay between Xt-k and xt. The figure 
is equivalent to the upper part of Fig. [1] the only difference 
being that blocks of size 7 are arranged as a two-dimensional 
grid. 




III. TWO-DIMENSIONAL REPRESENTATION 

The main purpose of the mapping described above is 
to arrange the reordered time series in a 1+1-dimensional 
plane in such a way that the interactions become local. 
Following Giacomelli and Politi [ll[ the time series is 
divided into equidistant segments of n elements which are 
arranged line by line on top of each other. This means 
that the index t is mapped to the position 

(:) = 0"f ) < s » 

Obviously this procedure turns the long delay n into a 
local interaction. For example, drawing the original time 
series x t in such a 1+1-dimensional representation the 
long delay n turns into a nearest-neighbor interaction in 
vertical direction while the short delay is still non-local 
(see Fig. [2]). However, drawing the reordered time se- 
ries y t in a 1+1-dimensional representation the long as 



FIG. 3: (Color online) Corresponding reordered time series in 
a two-dimensional representation according to Eq. (j4}, where 
the bulk interactions become local. 



well as the short delay become local, as shown in Fig. 03 
Note that the corresponding boundary conditions are not 
periodic but shifted in vertical direction, connecting sub- 
sequent blocks periodically in a spiral-like manner. 

Although the couplings in Fig. [3] are local, they are 
still biased towards north-east. Moreover, the coupling 
scheme exhibits vertical dislocation lines. As illustrated 
in Fig. 31 these irregularities can be removed by redraw- 
ing the figure in such a way that all updates have the 
same orientation in the x?/-plane. This representation 
allows one to relate the original time series with 1+1- 
dimensional cellular automata on a tilted square lattice. 
However, by rearranging the lattice one obtains skewed 
boundary conditions with a non-local delay connecting 
blocks with vertical distance n — k. 
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FIG. 4: (Color online) The same lattice as in Fig. [3] plotted 
in a way such that all updates in the bulk are geometrically 
identical. The boundary conditions one obtains can be seen 
as being time-delayed (time flows upwards). 

To demonstrate how the transformation works let us 
consider a simple deterministic update rule where x t G 
0, 1 is a binary time series. The update rule is given by 
the Boolean function 

Xt = Xt - k = {I 'iS;zt Xt ~ k • < 9 > 

where denotes a logical XOR. Such a binary time se- 
ries can be visualized by plotting x t as black and white 
pixels at position (x, y). The results are shown in Fig. [5l 
where we used the delays n = 300 and k = 227. Starting 
with a single non-zero entry in the initial state x n / 2 = 1 
the iteration of the update rule $1]) produces an irregular 
pattern of pixels, which is shown in the upper panel of 
Fig.0 Applying the transformation (|4]) the pixels are or- 
dered with a bias to the north-east, resulting in a tilted 
Sierpinsky gasket (see middle panel). Finally, plotting 
the same data in such a way that the tilt is removed 
one obtains the usual form of the Sierpinsky gasket (see 
lower panel). However, as exemplified by the arrow, the 
boundary conditions are no longer periodic, instead they 
involve a non-local shift 5 = n — k in vertical direction. 



IV. STOCHASTIC UPDATE RULE RELATED 
TO DIRECTED PERCOLATION 

Let us now turn to a simple but non-trivial example of 
a stochastic update rule. For a binary time series a prob- 
abilistic update (xt-ni%t-k) x t according to Eq. (pQ) 
is determined by 

f(l,Xi,Xj) = p XuXj (10) 

f(0jXi,Xj) = 1— p Xi ,Xj 

with four control paramters poo , poi , pio , pn . In what fol- 
lows let us assume that poo =0. In this case the time 
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FIG. 5: Visual representation of the time series generated 
by the update rule J9}. The three panels correspond to the 
examples shown above in Figs. [2] [3] and U (see text). 
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series consisting of zeroes is a fixed point of the dynam- 
ics. In non-equilibrium statistical physics such a config- 
uration which can be reached but not be left is called an 
absorbing state. Whether or not this absorbing state is 
stable against perturbations depends on the magnitude 
of the remaining control parameter poi,pio and p\\. For 
example, setting 

Poi=Pio=P, pu=2p-p 2 (11) 

and varying p between and 1 one observes the following 
phenomenological behavior: 

• If p is very small the time series quickly approaches 
the absorbing series consisting of zeroes. 

• For large p the dynamics approaches a fluctuating 
steady state with a non-vanshing stationary expec- 
tation value of x t . The probability to reach the 
absorbing configuration is very low and decreases 
with increasing n. 

• At a certain threshold p c ~ 0.6447 one observes a 
power-law decay of the density in a finite temporal 
range which grows with n. 

This behavior reminds of a non-equilibrium phase tran- 
sition from a fluctuating phase into an absorbing state. 
In fact, using the transformation (|2j) the update rule be- 
comes equivalent to that of a Domany-Kinzel cellular au- 
tomaton [l2j], which is known to exhibit a second-order 
phase transition belonging to the universality class of 
Directed Percolation (DP) 0, [3, EE O- Using the 
present notation the 1+1-dimensional Domany-Kinzel 
model is defined on a tilted square lattice with coor- 
dinates (x,y). Each lattice site can be either active 
(s(x,y) = 1) or inactive (s(x, y) =0). The model evolves 
by parallel updates, i.e. the new horizontal line at y + 1 
is obtained by setting 

Six V + 1) '= i 1 With P r0babiHt y Ps(x-l,y),s(x+l,y) 

^ ' 1 otherwise. 

(12) 



For the choice poi = Pio = P and pu = 2p — p 2 the 
Domany-Kinzel model reduces to directed bond percola- 
tion. This model is known to exhibit a continuous phase 
transition belonging to the universality class of directed 
percolation at the critical point p c = 0.6447001(2) if the 
system size is infinite. In fact, as shown in Fig. [6l the 
transformation (|4|) maps an apparently disordered time 
series into an ordered one, where typical DP cluster can 
be seen. 

It should be stressed that in the present model the 
corresponding DP process takes place on a finite lattice 
so that for any finite n there is no phase transition in 
a strict sense. Nevertheless it is possible to observe the 
typical signatures of DP critical behavior within a certain 
temporal range which grows with n, as will be shown in 
the following. 



A. Two-point correlation function 

In order to see a signature of DP critical behavior we 
tried to identify the critical exponents in 1+1 dimensions 

(3 = 0.276486(8), u± = 1.098654(4), v\\ = 1.733847(6). 

(13) 

To this end we iterated the time series slightly above the 
critical point p = p c + 0.001 = 0.6448 and measured the 
connected part of the two-point correlation function 

c(r) = (x t x t - T ) - (x t ) 2 . (14) 

Before the average was taken the time series was equili- 
brated over 2 x 10 9 iterations in order to reach a station- 
ary state. Moreover, we chose a very large delay n — 4096 
to prevent the system from entering the zero sequence 
due to finite-size effects. Figure [7] shows a train of spikes 
at regularly spaced times r s which are caused by corre- 
lations between subsequent rows in the two-dimensional 
representation. The asymptotic envelope of these spikes 
seems to obey an asymptotic power law 

c(r s ) ~ (15) 
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FIG. 7: Two-point correlation function of the time series c(r) 
denned in Eq. (|14[) . The dashed line demonstrates that the 

spike peaks decay as r s . The inset shows a zoom of the 
first and the last spike. 

with an exponent 7 = 0.32(4). The relation to DP pre- 
dicts this exponent to be given by 

7 = 2/3/i/y ~ 0.318. (16) 

A similar attempt to obtain the spatial correlation ex- 
ponent 2(3/ v± from the form of a single spike (see inset) 
fails. This can be explained as follows. In the central 
panel of Fig. [6] the form of the spike would correspond to 
a correlation function in horizontal direction, whereas in 
the symmetrized representation shown in the right panel 
this correlation function would be tilted. Therefore, the 
spike profile is given by an interplay of both exponents v\\ 
and z/j_, making it difficult to distinguish between them. 



B. Dynamical scaling of the pair connectedness 
function 

In order to identify the DP critical exponents more 
clearly, we measured the counterpart of the so-called pair 
connectedness function (see e.g. [14[) at criticality p = p c . 
The iteration starts with a single active seed x n = £ n , no . 
After reordering the time series and representing it in the 
right panel of Fig. [6l we measure the density p(A, y) in 
vertical distance y and horizontal distance A from the 
seed. Since the pair connectedness function is known to 
obey the scaling form 

P{^,y) = A" 2/3/iyj - h(y/A 1/z ) (17) 

with a universal scaling function h, the exponents can 
be determined by plotting p(A,y)A 2 P/ u± versus y/A 1 ^ 
in such a way that data sets for different values of A 
collapse. Plugging in the known exponents of DP one 
obtains a convincing data collapse, as shown in Fig. [HJ 



FIG. 8: (Color online) Data collapse for the counterpart of 
the pair connectedness function for n = 512 and k = 1. 

This confirms unambiguously the critical behavior we ex- 
pected to get. 



V. DISCUSSION 

In this paper we have shown how a time series evolv- 
ing according to a nonlocal update rule can be mapped 
onto a local process in two dimensions with special time- 
delayed boundary conditions. One interesting question 
is whether this result holds also for the continuous case. 
For example, discretizing a differential equation of the 
form 

x(t) = F[x(t)} + G[x(t - t)] (18) 

with delay r and arbitrary functions F and G one gets 

x m ~ Xm-1 h[F(x m ) + G(x m - k )} , (19) 

where h is the step size and k = r/h is the discrete ana- 
logue of the delay. In the present paper we have shown 
that the equation 

X m ~ Xm-s = h[F(x m ) + G(x m -k)} (20) 

with 1 < k < s exhibits (up to boundary conditions) 
the same properties as long as s and k are coprime. The 
question is whether it is possible to find an appropriate 
limit of (|2Q|) and recover ([T5|) . Work in this direction is 
currently under way [171 ]. 
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APPENDIX A: PROOF OF THE 
TRANSFORMATION 

The transformation (j4]) is proven in two steps. First 
we show that two sites of the time series separated by 
a time delay of At = n are still separated by the same 
delay after the transformation. Then we show that sites 
that were originally separated by a time delay At = k are 
mapped onto new sites either with a time delay At = 1 
or At = n + 1. 

We start with the sites separated by n time steps, that 
is t and t — n. According to the transformation rule we 
get 

T(t) -T(t-ri) = (tq) mod n + n[t\ n 

— (t — n)q mod n —n [(t — n)\ n 

= (tq) mod n 

= n([t] n - [(t-n)\ n ) = n, (Al) 



where we used the fact that for two sites of different 
blocks (which is always the case here) one has [t\ n — 
L(t-n)J n = l. 

Next we consider the case where sites have a delay of k, 
i.e. 

T(t) - T(t - k) = (tq) mod n - ((t - k)q) mod n 

+ n(\_t\ n - L*-lJn) (A2) 

Here we have to distinguish two cases. If t is a multiple 
of n this expression reduces to 

T(t)-T(t-k) = -(-kq) mod n+n([t\ n - [t - l\ n ) 

V v ' V v ' 

=n-l =1 

= 1 (A3) 
On the other hand, if t is not a multiple of n we get 



T(t) - T(t -k) = (tq) mod n - (tq - kq) mod n + n([t\ n - [t - lj n ) 

= (tq) mod n — {(tq) mod n — (kq) mod n^j mod n + n ( \t\ n — \ t — lj n ) 
e{i,...,n-i} e{i,...,n-i} =i 

' 1 

n - 



= 1 + Tl([t\ n - [t-l\ n ) = 



if 

1 if 



t mod n > k 
t mod n < k 



(A4) 
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